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We study the spontaneous formation of vortices during the superfluid condensation in a trapped 
fermionic gas subjected to a rapid thermal quench via evaporative cooling. Our work is based on the 
numerical solution of the time dependent crossover Ginzburg-Landau equation coupled to the heat 
diffusion equation. We quantify the evolution of condensate density and vortex length as a function 
of a crossover phase parameter from BCS to BEC. The more interesting phenomena occur somewhat 
nearer to the BEC regime and should be experimentally observable; during the propagation of the 
cold front, the increase in condensate density leads to the formation of supercurrents towards the 
center of the condensate as well as possible condensate volume oscillations. 

PACS numbers: 



Introduction. — Atomic gases undergoing a superfluid 
condensation are very clean and simple systems to ad- 
dress non-equilibrium quantum dynamics of supercon- 
ductors. As an example of important non-equilibrium 
behavior, the present paper focuses on the way in which 
a rapid quench from the normal phase leads to the for- 
mation of a random tangle of vortices. In addition to 
cosmology (lj, the character and dynamics of these de- 
fects is of great interest to a number of different physics 
subdisciplincs ranging from atomic phy sics, ^ |j| to con- 
densed matter || and high temperature superconduc- 
tivity To this end, theoretical investigations of 
rapid quenches of atomic Bose gases have been under- 
taken Q and appear consistent with experiments. In 
this Rapid Communication we make predictions for fu- 
ture observations on trapped Fermi gases throughout the 
entire crossover from BCS to BEC, with special emphasis 
on the unitary regime. 

These fermionic systems are in many ways more suit- 
able systems to study the dynamics of quantum gases. 
The Pauli principle suppresses three body collisions mak- 
ing it possible to study the strong interaction regime. For 
the Bose counterparts, collisions are frequent and lead 
to condensate decoherence. A related advantage of the 
Fermi gases is the opportunity to explore tunable interac- 
tion strengths which, in turn, affects the dynamics. This 
tunability is associated with the crossover between the 
BCS (where the inter- fermionic attraction is weak) and 
the Bosc-Einstcin condensed (BEC) limits (where the at- 
traction is strong). This crossover is entirely accessible 
through application of a magnetic field, and the exploita- 
tion of so-called Feshbach resonances 1,0. Our work is 
based on a numerical simulation of of the time-dependent 
Ginzburg-Landau equation (TDGLE) |ll|, [l^] in the pres- 
ence of thermal (white) noise x- This TDGLE equation, 
which in many ways is one of the most fundamental equa- 
tions in condensed matter physics, represents a differ- 



ential equation characterizing the dynamics and spatial 
dependence of the pairing gap parameter. Microscopi- 
cally, one can demonstrate that f§, @ the TDGLE 
coincides with the energy conserving Gross-Pitacvski de- 
scription in the strong attraction or BEC limit, where 
the dissipation is minimal Similarly in the weak attrac- 
tion limit the TDGLE leads to the well known diffusive 
dynamics of BCS theory. 

It should be stressed that TDGL theory refers to 
the gap or order parameter, and not to the underly- 
ing fermions. In this way the number of fermions, or 
the Fermi-Fermi interaction parameters do not enter, as 
they have been effectively integrated out. In contrast to 
static GL theory, which is a consequence of Gorkov the- 
ory, there is no fully rigorous derivation of this dynamics; 
alternative dynamical schemes such as time dependent 
Bogoliubov deGennes (TDBDG) theory [|l(| , which have 
similar antecedents in Gorkov theory, are also not rigor- 
ous. 

By contrast to this TDBDG approach, our studies are 
in the highly non-equilibrium regime. We stress that 
there are essentially no alternative effective simulation 
techniques for addressing the condensate evolution, asso- 
ciated with rapid quenches. Monte Carlo and first princi- 
ple quantum dynamics are limited to small system sizes, 
short times and typically zero temperature; moreover, 
they deal with the dynamics of single particles, not the 
collective dynamics of the order parameter (in the strong 
coupling limit) which we consider here. 

Theory. — We model the temperature equilibration 
in a physical way as a non-uniform process, associated 
with evaporative cooling. The effects of the BCS-BEC 
crossover enter most prominently via the change in the 
complex time relaxation rate of the TDGLE || [l3|, [b|. 
As a function of the continuous crossover from BCS to 
BEC we see that the number and lifetime of spontaneous 
vortices increases, so that the steady state condensate is 
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slower to form. Among the most interesting observations 
in this paper pertains to the near unitary regime (or mid- 
point between BCS and BEC) where the interactions are 
strongest and the TDGLE- based approach is the most 
appropriate. 

For notational purposes, it is convenient to introduce 
a parameter 9 = . . . n/2 that tunes the equation from 
BCS (6> = 0) to BEC (6 = tt/2). Physically this parame- 
ter can be viewed as representing an applied magnetic 
field in conjunction with a Feshbach resonance. This 
tunes the strength of interaction between the fundamen- 
tal fermionic particles (see e.g. ||, |l3|, [l4|] for a rigorous 
derivation) . 
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FIG. 1: (Color online) Isosurfaces and phase-cuts of the order 
parameter at different times. The upper row shows the 3D iso- 
surfaces for the BCS (6 = 0) case at times t — 0.8r, 0.9r, l.lr, 
where r is the time needed to reach half of the condensate 
steady state volume [see Fig. |^ (top)]. The phase-cuts (phase 
ifi = argi/>) are cuts through the center in the xy-plane - vor- 
tices are marked as the end-points of 27r-phase jump lines 
(sharp black- white transition). The lower row shows the same 
pictures, but towards the BEC limit represented by 6 = 85°. 

The TDGLE is given by 

e l() d t ip = A[r, T(r, t)]i> - |V>| 2 ^ + \ VV + x(r, t) (1) 

where ij) = ip{ r ,t) is the order parameter and \ is uni- 
formly distributed thermal noise with fluctuation tem- 
perature T x . 

The coefficient A[r, T(r, t)] takes into account the evap- 
orative cooling of the system and the trapping potential, 



i.e. A[r,T(r,t)) = [I - T{r,t)\ - U Q r 2 , where U is the 
trapping potential which is determined by the trap size. 
Eq. (|l|) is written in dimensionless units, with length scale 
£o , time scale To , and the order parameter scale ip® . Tem- 
peratures are measured in units of T c . For (fermionic) 
Li 6 these values are estimated as: £o = 3.2 • 10 _6 m, 
t = 1.1 • 10~ 2 s, and concentration ipQ = 2.9 • 10 16 m -3 . 
The dimensionless fluctuation temperature for a conden- 
sate at lOnK is in the interval T x G [10~ 4 ; 10~ 3 ], depend- 
ing on the condensed atoms. 

The trap is assumed to be spherically-symmetric. Thus 
we calculate the temperature profile T(r,t) due to evap- 
orative cooling by solving the spherically-symmetric heat 
diffusion equation 



d t T{r, t) = V 



(2) 



with an initial uniform temperature T(r < R, t = 0) = To 
and final temperature T{r > R,t) = Tf that is fixed at 
the boundary r = R. The heat diffusion constant V is 
rcnormalizcd by £o/ T o; typical values for Li 6 : T> ~ 3. 

All simulation runs start with random initial condi- 
tions in the normal phase (|^| 2 <C 1) at temperatures 
To larger than the critical temperature T c . The evapo- 
rative cooling is initiated at time t = on the surface 
of our spherical trap with radius 75£o- The cold front 
surface propagates towards the center and quenches the 
atomic gas below the critical temperature. The dimen- 
sionless diffusion constant was typically set to D = 10 
and the fluctuation strength T x = 1CP 4 . This mechanism 
of condensate nuclcation is to be contrasted with that de- 
scribed in Ref. where the cooling occurred uniformly 
in space, thus corresponding to an infinite diffusion con- 
stant. Another crucial difference with this previous work 
is the equilibration mechanism: the time evolution of the 
condensate used in was based on the Gross-Pitaevskii 
model where the dissipation was introduced by trunca- 
tion of high-frequency modes. By contrast, our TDGLE 
model includes a complex relaxation rate which is deter- 
mined by the crossover physics through the phase factor 
e 10 . Experiments on these cold gases [|| cannot probe 
very deeply into the BCS or BEC regimes, but are con- 
fined to the so-called "unitary" midpoint. Nevertheless, 
at unitarity, pairs are reasonably long lived so that 
we can associate 9 « 70° — 85° with the physically 
accessible regime. 

Numerical Results. — Our numerical calculations were 
done for volumes discretized in up to 512 3 grid points, 
averaged over up to 50 different initial conditions and 
time evolved for 65 different crossover phase values. The 
timestep in dimensionless units was chosen to be 0.1 and 
the total simulation time up to 3000 at ^-values close 
to the BEC limit. Our quasi-spectral split step method 
to solve the TDGLE which uses fast Fourier transforms, 
is much more stable than the traditional finite differ- 
ence method. As a result of employing modern graph- 
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ics processing units, our systems are an order of magni- 
tude larger than in recent work on the bosonic BEC [|| . 
This allows us to simulate more realistic physical situa- 
tions fn§. 

We begin with an illustration of the time evolved con- 
densation process for the BCS limit (9 = 0) and a near- 
BEC situation (9 = 85°). Throughout this paper we 
avoid the strict fermion-based BEC limit (9 = tt/2) since 
without dissipation the condensate does not form and nu- 
cleation of vortices is completely suppressed. By intro- 
ducing a complex relaxation rate in the TDGLE, we avoid 
having to include the interactions between condensed and 
non-condensed pairs which do not enter as naturally, in 
the fermionic- BEC limit. These were essential for ad- 
dressing the bosonic counterpart experiments ^ p7| . 

Plotted in Fig. [l] arc the isosurfaces for constant con- 
densate magnitude. We also show cross-sections in the 
xy-plane through the center of the system, indicating the 
phase of the order parameter ip = aigtp. These lat- 
ter plots appear below the counterpart isosurface pic- 
tures and contain information about spontaneous vor- 
tices. The isosurfaces for the condensate magnitude are 
color-coded according to fixed condensate density values 
p s = \ip\ 2 , where we used p s = 0.1,0.2,0.3. The three 
different panels from left to right correspond to three 
different times close to the typical time r at which the 
condensate forms, i.e. the time when the condensate vol- 
ume reaches half of its saturation value. 

One can see from Fig. [j], that most spontaneous vor- 
tices are connected to the surface of the condensate. 
Their time evolution is entirely accessible in our simu- 
lations Jl8| and one can follow their decay in the bulk 
and at the condensate surface. In the phase cross-section 
plots, these vortices are even more visible appearing as 
topological defects (phase singularities). Although vor- 
tices are most plentiful at the surface of the condensate, 
we find their decay occurs mostly in the bulk. 

It is clear that in the BCS limit the condensate forms 
quickly occupying the entire trap volume in a relatively 
short time period. At much longer time scales (than 
shown here), the condensate will appear uniform. The 
near-BEC limit exhibits an interesting contrast: Here the 
condensate expands slowly, and vortices are much more 
plentiful, decaying even more gradually. An interesting 
feature is shown in the middle panel of the last row of 
Fig. [l| We observe in the course of condensate formation, 
concentric rings in the phase cross-section interrupted by 
trapped topological defects. This corresponds to the gen- 
eration of supercurrents from the surface of the conden- 
sate towards its center. The phenomenon of transient su- 
percurrent J s = p s Vtp generation can be understood as 
follows: In the BEC limit the TDGLE can be written in 
the hydrodynamic form ]l9| with the condensate density 
p s satisfying the continuity equation dp s /dt = — VJ S . 
During the quench, p s changes from zero to its maximal 
value leading to the formation of supercurrents. 
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FIG. 2: (Color online) a) Condensate volume V c /Voo normal- 
ized to the asymptotic condensate volume V x with simulation 
time t/ro for different crossover phases 9. The horizontal line 
at Vc/Voo = 1/2 can be used to read the typical condensate 
formation time r. Inset: the behavior of the condensate vol- 
ume near the "overshoot" for 9 — 85° for the full simulation 
(labeled 3D) and the spherically-symmetric equation (labeled 
ID), b) Time evolution of the vortex length L v /Voo for differ- 
ent 9. Inset: a log-log representation illustrating a power-law 
decay for different 9. The straight dashed line corresponds to 
a 1/t decay. 



Figure |2j presents a plot of the time evolution of the 
condensate volume V c , as well as the vortex length L v 
for different values of the crossover phase 9. Both quan- 
tities are normalized to the asymptotic final condensate 
volume. Here the condensate volume is calculated from 
the number of grid points with condensate density larger 
than 0.2. Assuming a spherical shape of the condensate 
(reflecting the trapping potential) we establish the diam- 
eter of the condensate; we then calculate the number of 
grid points with p s < 0.2 inside this sphere to estimate 
the vortex length. In the BCS limit the dissipation is 
maximum which is responsible for the rapid formation of 
the condensate. Increasing 9 as one approaches the near- 
BEC gradually delays the formation of the condensate 
leading to increase of the maximum total vortex length. 

The vortices show a roughly 1 /i-decay of the total vor- 
tex length jl], 0], In Fig. |^b the time dependence of the 
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total vortex length is plotted from the BCS to the near- 
BEC limit. The inset shows some of these same curves on 
a log-log scale and compares with a linear curve (dashed) 
representing a 1/t decay. In two dimensions, bulk vor- 
tices annihilate in pairs, the relaxation should be propor- 
tional to the square of the vortex concentration, i.e. the 
relaxation follows a power-law pc| . Similarly, a 1/t is 
expected for the decay of vortex lines in 3D for the BCS 
case and 1/i 3 / 2 near the BEC limit [21 . The deviations 
from 1/t power law are likely caused by finite size effects. 

The condensate volume shows a well pronounced over- 
shoot with subsequent oscillations in time. This appears 
around unitarity, albeit closer to the BEC limit. This 
overshoot is not generic, since it disappears in the BCS 
limit, presumably as a consequence of the large dissi- 
pation. In the near-BEC limit, where the dissipation 
is strongly suppressed, the condensate takes longer to 
form than the time associated with a thermal quench. 
We quantify this "condensate breathing" by studying the 
height of this condensate volume overshoot in the inset 
of Fig. H Plotted here is the maximum of the volume of 
the condensate as a function of 9. Also indicated is the 
maximum vortex length within the condensate. It can be 
seen that the maximum of the condensate overshoot oc- 
curs at 6 = 82.5°. The appearance of the overshoot is also 
characterized by the time scale t c ,max, where the conden- 
sate reaches its maximum, which diverges below 9 ~ 60° 
(since there is no overshoot) and near the BEC limit (be- 
cause the time of the occurrence of the overshoot diverges 
for 9 — > 7r/2). A minimum of i c ,max is observed at 75°, 
i.e., the overshoot happens fastest. Fig. || indicates the 
typical time scale of the condensate formation r and the 
nearly identical time scale t v ,max, which is the time at 
which the vortex length is maximum. Both these times 
generally depend only weakly on 9, except near 6 — > tt/2. 
Also the typical decay time of the vortices t v . decay ap- 
pears to be nearly independent of 8 below ~ 85° (see 
Fig. ||b). Using the experimental time scales for Li 6 , we 
can conclude from Figs. ||& [| that the condensate forms 
within a few seconds and vortices annihilate on the same 
time scale. 

The overshoot and subsequent time oscillations of the 
condensate may be related to excitation of sound waves 
produced in a spherical trap by the quench. To fur- 
ther clarify this phenomenon we solved the spherically- 
symmetric TDGLE and obtained similar behavior as 
shown in the inset to Fig. ||a; the main difference being 
a delay in the condensate nucleation |2^] . In the process 
we verified that the period of these oscillations is close 
to the time needed for sound waves to traverse the con- 
densate. Indeed in our units the speed of sound (which 
emerges from the BEC limit of the TDGLE) is equal to 
unity in the long-wavelength limit, while the oscillation 
period obtained from Fig. ||a is approximately 140 time 
units. This is close to the time needed for the sound wave 
to cross our system with diameter 150 length units. Us- 
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FIG. 3: (Color online) Time scales extracted from the time 
evolution of the condensate volume and vortex length [time 
in units of To], t is the time needed to reach half of the 
steady state condensate volume. f c , max and i„ imax ar e the 
times when the condensate volume and vortex length has a 
local maximum respectively, i^decay is the time over which 
the vortex length decays to 1/e of it maximum value. Inset: 
Maximum of the volume of the condensate V c ,max and vortex 
length L VtTna x normalized to Voo versus the crossover phase 9. 
The data is thermally averaged over 10 to 50 realizations. 



ing the value for £o and To for Li 6 , the speed of the second 
sound is on the order 10~ 4 m/s, comparable to the mean 
velocity of the atoms at lOnK. 

There has been recent excitement over numerical ap- 
proaches for addresses vortices in Fcrmionic supcrfluids, 
particularly with the implementation [[l6| of a time de- 
pendent Bogoliubov dcGcnnes numerical scheme. Unlike 
TDGLE, this approach addresses the fermionic degrees 
of freedom in conjunction with the fcrmionic gap parame- 
ter. Generally it has been applied to T = 0. Importantly, 
it does not incorporate the fact that the superconducting 
order parameter need not be the same as fermionic exci- 
tation gap. Related to this last observation is the fact 
that these schemes do not accomodate non-condensed 
bosons which are expected to enter into the dynamics 
away from the strict BCS regime. As a consequence, 
our understanding of vortices in fermionic superfluids will 
require the exploration of a number of different numeri- 
cal approaches, including the TDGLE scheme we address 
here. 
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